#install.packages("rgbif")
#install.packages("ggplot2")
#install.packages("raster")
#install.packages("rasterVis")
#install.packages("maps")
#install.packages("dismo")
#install.packages("rgdal")
library(rgbif)
Species.Name <- "Vulpes macrotis"
#Search GBIF and download observations
Species.key <- name_backbone(name=Species.Name)$speciesKey
locations <- occ_search(taxonKey=Species.key,
decimalLatitude = '30,50',
decimalLongitude='-130,-75',
hasCoordinate = TRUE,
return='data', limit=5000)
#Remove localities with uncertain coordinates - meaning the observer entered approximate latitude and longitude data. This is common with very old specimens/records that were collected prior to handheld GPS devices being invented
locations <- subset(locations, coordinateUncertaintyInMeters <= 5000)
head(locations)
## # A tibble: 6 x 148
## key scientificName decimalLatitude decimalLongitude issues datasetKey
## <chr> <chr> <dbl> <dbl> <chr> <chr>
## 1 2557… Vulpes macrot… 33.2 -116. cdrou… 50c9509d-…
## 2 2557… Vulpes macrot… 33.2 -116. cdrou… 50c9509d-…
## 3 2563… Vulpes macrot… 33.0 -115. cdrou… 50c9509d-…
## 4 2574… Vulpes macrot… 32.5 -104. cdrou… 50c9509d-…
## 5 2597… Vulpes macrot… 33.0 -115. cdrou… 50c9509d-…
## 6 2005… Vulpes macrot… 33.3 -116. cdrou… 50c9509d-…
## # … with 142 more variables: publishingOrgKey <chr>, installationKey <chr>,
## # publishingCountry <chr>, protocol <chr>, lastCrawled <chr>,
## # lastParsed <chr>, crawlId <int>, extensions <chr>, basisOfRecord <chr>,
## # taxonKey <int>, kingdomKey <int>, phylumKey <int>, classKey <int>,
## # orderKey <int>, familyKey <int>, genusKey <int>, speciesKey <int>,
## # acceptedTaxonKey <int>, acceptedScientificName <chr>, kingdom <chr>,
## # phylum <chr>, order <chr>, family <chr>, genus <chr>, species <chr>,
## # genericName <chr>, specificEpithet <chr>, taxonRank <chr>,
## # taxonomicStatus <chr>, dateIdentified <chr>,
## # coordinateUncertaintyInMeters <dbl>, stateProvince <chr>, year <int>,
## # month <int>, day <int>, eventDate <chr>, modified <chr>,
## # lastInterpreted <chr>, references <chr>, license <chr>, identifiers <chr>,
## # facts <chr>, relations <chr>, geodeticDatum <chr>, class <chr>,
## # countryCode <chr>, recordedByIDs <chr>, identifiedByIDs <chr>,
## # country <chr>, rightsHolder <chr>, identifier <chr>,
## # http...unknown.org.nick <chr>, verbatimEventDate <chr>, datasetName <chr>,
## # gbifID <chr>, collectionCode <chr>, verbatimLocality <chr>,
## # occurrenceID <chr>, taxonID <chr>, catalogNumber <chr>, recordedBy <chr>,
## # http...unknown.org.occurrenceDetails <chr>, institutionCode <chr>,
## # rights <chr>, eventTime <chr>, occurrenceRemarks <chr>,
## # identificationID <chr>, name <chr>, infraspecificEpithet <chr>,
## # informationWithheld <chr>, identificationRemarks <chr>,
## # identifiedByIDs.type <chr>, identifiedByIDs.value <chr>,
## # X.1f2c0cbe.40df.43f6.ba07.e76133e78c31. <chr>, individualCount <int>,
## # sex <chr>, continent <chr>, institutionID <chr>, dynamicProperties <chr>,
## # county <chr>, identificationVerificationStatus <chr>, language <chr>,
## # type <chr>, locationAccordingTo <chr>, preparations <chr>,
## # identifiedBy <chr>, georeferencedDate <chr>, higherGeography <chr>,
## # nomenclaturalCode <chr>, georeferencedBy <chr>, georeferenceProtocol <chr>,
## # georeferenceVerificationStatus <chr>, locality <chr>,
## # verbatimCoordinateSystem <chr>, otherCatalogNumbers <chr>,
## # organismID <chr>, previousIdentifications <chr>,
## # identificationQualifier <chr>, accessRights <chr>, collectionID <chr>, …
locations.sub <- locations[,c("species", "decimalLongitude", "decimalLatitude")]
#rename the headers on each column to make our commands easier below
colnames(locations.sub) <- c("Species", "Longitude", "Latitude")
#Save dataset as a .csv file on your computer
write.csv(locations.sub, "Locations_KF.csv", quote=F, row.names=F)
library(ggplot2)
# Download Country Boundaries
#install.packages("maps")
world <- map_data("world")
map <- ggplot()+
geom_polygon(data=world, aes(x=long, y=lat, group=group), color='grey40', fill='transparent', size=0.25)+
geom_point(data=locations.sub, aes(x=Longitude, y=Latitude), size = 2, alpha =0.75) +
theme_bw()+
coord_fixed(xlim = c(-125, -110), ylim = c(32, 43), ratio = 1.3)
map
ggsave("LocationsMap.jpg", height=5, width=7, units="in", dpi=600)
library(raster)
## Loading required package: sp
#describe the extent for which we want to download data
e <- extent(c(-125, -110, 32, 43)) ## min longitude, max longitude, min latitude, max latitude
#download the data for our region of interest
bioclim <- getData('worldclim', var='bio', res=2.5) ## change this command from what we did in lab 8, we have increased the resolution, which allows us to download the full global dataset
#crop the dataset to just the coordinates we are interested in
predictors <- crop(bioclim, e)
#plot all of the variables
plot(predictors)
#create a folder called "Predictors"
dir.create("Predictors", showWarnings = F)
#run a for loop to save all 19 .asc files
for(i in 1:19){
file.name <- paste0("Predictors/", names(predictors[[i]]), ".asc")
writeRaster(predictors[[i]],file.name, overwrite=T)
}
dir.create("Maxent Output", showWarnings = F)
library(rasterVis)
## Loading required package: lattice
## Loading required package: latticeExtra
##
## Attaching package: 'latticeExtra'
## The following object is masked from 'package:ggplot2':
##
## layer
#Change your Species.Name object so that it has an underscore rather than a space in between genus and species (E.g, Canis latrans changes to Canis_latrans). This will make it possible to opben the .asc file that Maxent created
Species_Name <- gsub(" ", "_", Species.Name)
SDM.1 <- raster(paste0("Maxent Output/",Species_Name, ".asc"))
##Plot your SDM on a map with the outlines to the countries and the observation localities
map <- gplot(SDM.1, maxpixels = 100000) +
geom_tile(aes(fill = value))+
scale_fill_viridis_c(option = "C", na.value = "transparent") +
theme_bw()+
coord_fixed(xlim = c(-125, -110), ylim = c(32, 43), ratio = 1.3)+
xlab("Longitude")+
ylab("Latitude")
map
ggsave("SDM_KitFox_Current.jpg", height=5, width=7, units="in", dpi=600)
#install.packages("rgdal")
library(raster)
predictors.26 <- getData('CMIP5', var='bio', res=2.5, rcp=26, model='MI', year=50)
predictors.26 <- crop(predictors.26, e)
names(predictors.26) <- names(predictors)
plot(predictors.26)
### Save RCP 2.6 Predictor variables
#create a folder called "Predictors"
dir.create("Data Layers RCP 26", showWarnings = F)
#run a for loop to save all 19 .asc files
for(i in 1:19){
file.name <- paste0("Data Layers RCP 26/", names(predictors.26[[i]]), ".asc")
writeRaster(predictors.26[[i]],file.name, overwrite=T)
}
predictors.85 <- getData('CMIP5', var='bio', res=2.5, rcp=85, model='MI', year=50)
predictors.85 <- crop(predictors.85, e)
names(predictors.85) <- names(predictors)
plot(predictors.85)
### Save RCP 8.5 Predictor variables
#create a folder called "Predictors"
dir.create("Data Layers RCP 85", showWarnings = F)
#run a for loop to save all 19 .asc files
for(i in 1:19){
file.name <- paste0("Data Layers RCP 85/", names(predictors.85[[i]]), ".asc")
writeRaster(predictors.85[[i]],file.name, overwrite=T)
}
library(rasterVis)
library(ggplot2)
Species_Name <- gsub(" ", "_", Species.Name)
SDM.26 <- raster(paste0("Maxent Output RCP 26/",Species_Name, "_Data Layers RCP 26.asc"))
##Plot SDM on a map with the outlines to the countries and the observation localities
map <- gplot(SDM.26, maxpixels = 100000) +
geom_tile(aes(fill = value))+
scale_fill_viridis_c(option = "C", na.value = "transparent") +
theme_bw()+
coord_fixed(xlim = c(-125, -110), ylim = c(32, 43), ratio = 1.3)+
xlab("Longitude")+
ylab("Latitude")
map
ggsave("SDM_RCP26_KitFox.jpg", height=5, width=7, units="in", dpi=600)
library(rasterVis)
Species_Name <- gsub(" ", "_", Species.Name)
SDM.85 <- raster(paste0("Maxent Output RCP 85/",Species_Name, "_Data Layers RCP 85.asc"))
##Plot SDM on a map with the outlines to the countries and the observation localities
map <- gplot(SDM.85, maxpixels = 100000) +
geom_tile(aes(fill = value))+
scale_fill_viridis_c(option = "C", na.value = "transparent") +
theme_bw()+
coord_fixed(xlim = c(-125, -110), ylim = c(32, 43), ratio = 1.3)+
xlab("Longitude")+
ylab("Latitude")
map
ggsave("SDM_RCP85_KitFox.jpg", height=5, width=7, units="in", dpi=600)
#Compare niche overlap between current SDM and future SDMs under RCP conditions as well as future SDMs to each other for kit fox populations
#install.packages("dismo")
library(dismo)
#Calculate niche overlap using Schoener's D.
nicheOverlap(SDM.1, SDM.26, stat = "D")
## [1] 0.6881662
nicheOverlap(SDM.1, SDM.85, stat = "D")
## [1] 0.6535981
nicheOverlap(SDM.26, SDM.85, stat = "D")
## [1] 0.8335757
#Values close to 1 = a lot of overlap, values close to 0 = little overlap
Species.Name <- "Vulpes vulpes"
library(rgbif)
#Search GBIF and download observations
Species.key <- name_backbone(name=Species.Name)$speciesKey
locations <- occ_search(taxonKey=Species.key,
decimalLatitude = '32,43',
decimalLongitude='-125,-110',
hasCoordinate = TRUE,
return='data', limit=5000)
#Remove localities with uncertain coordinates
locations <- subset(locations, coordinateUncertaintyInMeters <= 5000)
head(locations)
## # A tibble: 6 x 145
## key scientificName decimalLatitude decimalLongitude issues datasetKey
## <chr> <chr> <dbl> <dbl> <chr> <chr>
## 1 2573… Vulpes vulpes… 38.4 -123. cdrou… 50c9509d-…
## 2 2596… Vulpes vulpes… 41.1 -112. cdrou… 50c9509d-…
## 3 2574… Vulpes vulpes… 40.7 -112. cdrou… 50c9509d-…
## 4 2574… Vulpes vulpes… 40.7 -112. cdrou… 50c9509d-…
## 5 2574… Vulpes vulpes… 40.7 -112. cdrou… 50c9509d-…
## 6 2576… Vulpes vulpes… 38.5 -122. gass84 50c9509d-…
## # … with 139 more variables: publishingOrgKey <chr>, installationKey <chr>,
## # publishingCountry <chr>, protocol <chr>, lastCrawled <chr>,
## # lastParsed <chr>, crawlId <int>, extensions <chr>, basisOfRecord <chr>,
## # taxonKey <int>, kingdomKey <int>, phylumKey <int>, classKey <int>,
## # orderKey <int>, familyKey <int>, genusKey <int>, speciesKey <int>,
## # acceptedTaxonKey <int>, acceptedScientificName <chr>, kingdom <chr>,
## # phylum <chr>, order <chr>, family <chr>, genus <chr>, species <chr>,
## # genericName <chr>, specificEpithet <chr>, taxonRank <chr>,
## # taxonomicStatus <chr>, dateIdentified <chr>, stateProvince <chr>,
## # year <int>, month <int>, day <int>, eventDate <chr>, modified <chr>,
## # lastInterpreted <chr>, references <chr>, license <chr>, identifiers <chr>,
## # facts <chr>, relations <chr>, geodeticDatum <chr>, class <chr>,
## # countryCode <chr>, recordedByIDs <chr>, identifiedByIDs <chr>,
## # country <chr>, rightsHolder <chr>, identifier <chr>,
## # http...unknown.org.nick <chr>, verbatimEventDate <chr>, datasetName <chr>,
## # gbifID <chr>, collectionCode <chr>, verbatimLocality <chr>,
## # occurrenceID <chr>, taxonID <chr>, catalogNumber <chr>, recordedBy <chr>,
## # http...unknown.org.occurrenceDetails <chr>, institutionCode <chr>,
## # rights <chr>, eventTime <chr>, identificationID <chr>, name <chr>,
## # coordinateUncertaintyInMeters <dbl>, occurrenceRemarks <chr>,
## # informationWithheld <chr>, recordedByIDs.type <chr>,
## # recordedByIDs.value <chr>, identifiedByIDs.type <chr>,
## # identifiedByIDs.value <chr>, infraspecificEpithet <chr>, sex <chr>,
## # continent <chr>, dynamicProperties <chr>, county <chr>, language <chr>,
## # type <chr>, preparations <chr>, reproductiveCondition <chr>,
## # recordNumber <chr>, georeferencedDate <chr>, higherGeography <chr>,
## # georeferencedBy <chr>, georeferenceProtocol <chr>,
## # georeferenceVerificationStatus <chr>, locality <chr>, startDayOfYear <chr>,
## # higherClassification <chr>, georeferenceSources <chr>,
## # X.1f2c0cbe.40df.43f6.ba07.e76133e78c31. <chr>, individualCount <int>,
## # elevation <dbl>, elevationAccuracy <dbl>, institutionID <chr>,
## # identificationVerificationStatus <chr>, locationAccordingTo <chr>,
## # identificationRemarks <chr>, …
locations.sub <- locations[,c("species", "decimalLongitude", "decimalLatitude")]
#Rename the headers on each column to make our commands easier below
colnames(locations.sub) <- c("Species", "Longitude", "Latitude")
#Save dataset as a .csv file on our computers
write.csv(locations.sub, "Locations_SN.csv", quote=F, row.names=F)
library(ggplot2)
map <- ggplot()+
geom_polygon(data=world, aes(x=long, y=lat, group=group), color='grey40', fill='transparent', size=0.25)+
geom_point(data=locations.sub, aes(x=Longitude, y=Latitude), size = 2, alpha =0.75) +
theme_bw()+
coord_fixed(xlim = c(-125, -110), ylim = c(32, 43), ratio = 1.3)
map
ggsave("Locations Map for Sierra Nevada Fox.jpg", height=5, width=7, units="in", dpi=600)
dir.create("Maxent Output 2", showWarnings = F)
library(rasterVis)
#Change our Species.Name object so that it has an underscore rather than a space in between genus and species. This will make it possible to opben the .asc file that Maxent created
Species_Name <- gsub(" ", "_", Species.Name)
SDM <- raster(paste0("Maxent Output 2/",Species_Name, ".asc"))
##Plot our SDM on a map with the outlines to the countries and the observation localities
map <- gplot(SDM, maxpixels = 100000) +
geom_tile(aes(fill = value))+
scale_fill_viridis_c(option = "C", na.value = "transparent") +
theme_bw()+
coord_fixed(xlim = c(-125, -110), ylim = c(32, 43), ratio = 1.3)+
xlab("Longitude")+
ylab("Latitude")
map
ggsave("Sierra Nevada Red Fox.jpg", height=5, width=7, units="in", dpi=600)
##Plot SDM26 on a map with the outlines to the countries and the observation localities
SDM26 <- raster(paste0("Maxent Output RCP 26 2/",Species_Name, "_Data Layers RCP 26.asc"))
map <- gplot(SDM26, maxpixels = 100000) +
geom_tile(aes(fill = value))+
scale_fill_viridis_c(option = "C", na.value = "transparent") +
theme_bw()+
coord_fixed(xlim = c(-125, -110), ylim = c(32, 43), ratio = 1.3)+
xlab("Longitude")+
ylab("Latitude")
map
ggsave("SDM26_SierraFox.jpg", height=5, width=7, units="in", dpi=600)
##Plot SDM85 on a map with the outlines to the countries and the observation localities
SDM85 <- raster(paste0("Maxent Output RCP 85 2/",Species_Name, "_Data Layers RCP 85.asc"))
map <- gplot(SDM85, maxpixels = 100000) +
geom_tile(aes(fill = value))+
scale_fill_viridis_c(option = "C", na.value = "transparent") +
theme_bw()+
coord_fixed(xlim = c(-125, -110), ylim = c(32, 43), ratio = 1.3)+
xlab("Longitude")+
ylab("Latitude")
map
ggsave("SDM85_SierraFox.jpg", height=5, width=7, units="in", dpi=600)
library(dismo)
#calculate niche overlap using Schoener's D.
nicheOverlap(SDM, SDM26, stat = "D")
## [1] 0.5812436
nicheOverlap(SDM, SDM85, stat = "D")
## [1] 0.5428118
nicheOverlap(SDM26, SDM85, stat = "D")
## [1] 0.8277318
#values close to 1 = a lot of overlap, values close to 0 = little overlap